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Abstract 

Least  squares  methods  based  on  first-order  systems  have  been  recently 
proposed  and  analyzed  for  second-order  elliptic  equations  and  systems.  They 
produce  symmetric  and  positive  definite  discrete  systems  by  using  standard 
finite  element  spaces,  which  are  not  required  to  satisfy  the  inf-sup  condition. 
In  this  paper,  several  domain  decomposition  algorithms  for  these  first-order 
least  squares  methods  are  studied.  Some  representative  overlapping  and 
substructuring  algorithms  are  considered  in  their  additive  and  multiplicative 
variants.  The  theoretical  and  numerical  results  obtained  show  that  the  clas¬ 
sical  convergence  bounds  (on  the  iteration  operator)  for  standard  Galerkin 
discretizations  are  also  valid  for  least  squares  methods. 
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1  Introduction 


Least  squares  methods  have  been  proposed  in  recent  years  for  second-order  el¬ 
liptic  problems,  Stokes  and  Navier-Stokes  equations;  see  Chang  [10],  Bochev 
and  Gunzburger  [2],  Pehlivanov,  Carey,  and  Lazarov  [15],  Cai,  Lazarov,  Man- 
teuflel,  and  McCormick  [5],  Cai,  Manteuffel,  and  McCormick  [7],  Bramble, 
Lazarov,  and  Pasciak  [3],  Bramble  and  Pasciak  [4],  Carey,  Pehlivanov,  and 
Vassilevski  [8],  Cai,  Manteuffel,  and  McCormick  [6],  Bochev,  Cai,  Manteuffel, 
and  McCormick  [1],  and  the  references  therein. 

Among  the  possible  approaches,  we  follow  here  the  one  introduced  in  the 
very  recent  works  of  Pehlivanov,  Carey,  and  Lazarov  [15]  and  Cai,  Manteuf¬ 
fel,  and  McCormick  [7].  The  second-order  elliptic  problem  is  rewritten  as  a 
first-order  system  and  a  least  squares  functional  is  introduced.  The  result¬ 
ing  discrete  minimization  problem  is  associated  with  a  bilinear  form  which  is 
continuous  and  elliptic  on  an  appropriate  space.  Therefore,  the  inf-sup  condi¬ 
tion  is  avoided  and  standard  finite  element  spaces  can  be  used.  The  resulting 
linear  system  is  symmetric,  positive  definite  and  has  condition  number  of  the 
same  order  as  standard  Galerkin  finite  element  stiffness  matrices,  0(l/fi^). 
An  interesting  alternative  approach  by  Bramble,  Lazarov,  and  Pasciak  [3]  is 
based  on  replacing  one  of  the  L^-terms  in  the  least  squares  functional  by  a 
discrete  L7“^-norm.  We  will  not  consider  here  such  an  alternative. 

The  goal  of  this  paper  is  to  extend  to  these  least  squares  methods  some  of 
the  classical  domain  decomposition  algorithms  which  have  been  successfully 
employed  for  standard  Galerkin  finite  elements.  We  show  that  optimal  and 
quasi-optimal  convergence  bounds  follow  easily  from  the  standard  Galerkin 
case.  Therefore,  domain  decomposition  provides  highly  parallel  and  scalable 
solvers  also  for  first-order  system  least  square  discretizations.  An  overview  of 
domain  decomposition  methods  can  be  found  in  the  review  papers  by  Chan 
and  Mathew  [9],  Dryja,  Smith,  and  Widlund  [11],  Dryja  and  Widlund  [13], 
and  Le  Tallec  [14]. 

This  paper  is  organized  as  follows.  In  the  next  section,  we  briefly  review 
the  first-order  system  least  squares  methodology  and  the  main  results  from 
[7].  In  Section  3,  we  introduce  and  analyze  our  domain  decomposition  algo¬ 
rithms;  overlapping  additive  Schwarz  methods  (with  coupled  and  uncoupled 
subspaces;  see  3.1),  overlapping  multiplicative  Schwarz  methods  (3.2),  and 
an  iterative  substructuring  method  (3.3).  In  Section  4,  we  present  numerical 
results  in  the  plane  that  confirm  the  theoretical  bounds  obtained. 
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2  First-Order  System  Least  Squares 


We  consider  the  following  second-order  elliptic  problem  on  a  bounded  domain 
ncR'^  ox 

'  -V  •  (AVp)  +  Xp  =  f  in 

<  p  =  0  on  Td,  (1) 

n-AVp  =  0  on 

Here  A  is  a  symmetric  and  uniformly  positive  definite  matrix  with  entries  in 
X  is  a  first-order  linear  differential  operator,  Fd  U  Fat  =  and  n 
is  the  outward  normal  unit  vector  to 

Defining  the  new  flux  variable  u  =  —  AVp,  the  system  (1)  can  be  rewritten 
as  a  first-order  system: 


u  +  AVp  =  0  in  D, 

V  •  u  -I-  Xp  —  f  in  fl, 

p  =  0  on  Td, 

n  •  u  =  0  on  Fat, 

This  system  can  be  extended  to  the  equivalent  system 


u  -|-  AVp  =  0 

in  fl. 

V-u-l-Ap  =  / 

in  fl. 

V  X  =  0 

in  fl. 

p  =  0 

on  Fd, 

n  •  u  =  0 

on  Tn, 

7^(A~^u)  =  0 

on  Fa), 

(2) 


(3) 


where  V  x  is  the  curl  operator  (in  two  dimensions  Vxu  =  0  means  = 

0)  and  7tU  =  u  x  n  (in  two  dimensions  7tU  =  u  •  r). 

The  following  least  squares  functionals,  Go  for  the  system  (2)  and  G  for 
the  augmented  system  (3),  were  studied  in  [5]  ([15]  for  the  case  A  =  0)  and 
[7],  respectively: 

Go(v,  p;  /)  =  llv  -f  AVg||i2(n)  -1-  ||V  •  v  -|-  -  /|lA2(fi)  (4) 


V(v,  9)  €  Wo(diu;fl)  X  y,  and 

G{'v,q;f)  =  ||v  +  AVg||i2(n)  +  ||V-v  +  A?-/l|A2(n)4-||Vx(A  V)||^2(q)  (5) 
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V(v,9)e  W  X  K 

Here  the  functional  spaces  are  defined  as 


Wo(dfu;0)  =  {v  G  H(diu;  fi)  :  n  •  v  =  0  on  Fjv}, 

Wo(curM;0)  =  {v  e  H(cur/yl;0)  :  7t(A“^v)  =  0  on  To}, 

W  =  Wo{div,  fi)  n  Wo(cMr/A;  fi), 

V  =  {q  e  :  ^  =  0  on  Fc}. 

The  least  squares  minimization  problems  for  (2)  and  (3)  are,  respectively: 
Find  (u,p)  €  Wo{div,Q,)  x  V  such  that 


Go(u,p;  /) 


inf 

{v,q)e'Wo{div;Q)xV 


Go(v,g;/); 


Find  (u,p)  €  W  xV  such  that 


(6) 


Simple  calculations  show  that  the  associated  variational  problems  are, 
respectively: 

Find  (u,p)  G  Wo(dfi;;  fi)  x  F  such  that 

«o(u,  p;  V,  q)  =  F{v,  q)  V(v,  q)  G  Wo(df u;  0)  X  F;  (8) 

Find  (u,p)  €  W  x  F  such  that 

«(u,p;v,q)  =  ^(v,^)  V(v,q)GWxF.  (9) 

Here  the  bilinear  forms  are 


Go(u,p;  v,§)  =  (u  +  AVp,v  +  +  (V  •  u  +  Xp,  V  •  \  -\-Xq)i,2, 

o(u,p;v,9)  =  ao(u,p;v,5)  +  (V  x  (A“^u),V  x  {A~W))l2 
and  the  right-hand  side  is 

-P'(v,g)  =  (/,V-v-l-Xg)j,2. 

In  [5],  it  was  proved  that  ao{v,q-,v,q)  is  equivalent  to  (continuous  and 
elliptic  with  respect  to)  the  'H{div;  H)  xH^{U)  norm  on  Wo{div,  Cl)xV,  under 
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the  assumption  that  a  Poincare- Friedrichs  inequality  holds  for  p  (denoted  by 
assumption  AO).  For  the  case  X  =  0,  this  was  already  proved  in  [15]. 

In  [7],  it  was  proved  that  a{v,q-,v,q)  is  equivalent  to  the  [li{div;Q.)  n 
H(cur/A;n)]  x  Ar^(O)  norm  on  W  x  V,  under  the  same  assumption  AO. 
Moreover,  under  three  additional  technical  assumptions  denoted  by  Al,  A2, 
A3,  it  is  proven  in  [7]  that  a(v,5;v,g)  is  equivalent  to  the  norm 

on  W  X  V  (d  =  2  or  3): 

Theorem  1  Let  b{u,p;v,q)  =  (u,  v)/jri(f2)d  -t-  (p,q)H'^(Q)  be  the  bilinear  form 
associated  with  the  norm. 

If  the  assumptions  A0-A3  of  [7]  are  verified,  then  there  exist  positive 
constants  a  and  /d  such  that 

a6(v,  q;  v,  q)  <  a(v,  g;  v,  q)  V(v,  g)  €  W  X  V, 

and 

a(u,p;v,5)  < 

V(u,p),(v,g)  e  w  X  y. 

Because  of  the  equivalence  of  a(-,  •)  and  b{-,-),  from  now  on  we  will  concen¬ 
trate  on  the  variational  problem  (9)  associated  with  the  augmented  system 

(3). 

We  introduce  a  triangulation  th  of  and  associated  finite  element  sub¬ 
spaces  X  C  W  X  V.  We  then  obtain  a  finite  element  discretization 

of  (9): 

Find  {uh,Ph)  ^  x  V  such  that 

a{uh,Ph',Vhjqh)  =  Fi'^h.qh)  V(v/„ g/i)  e  x  V'*.  (10) 

For  simplicity,  we  consider  continuous  piecewise  linear  finite  elements: 

=  {v  e  CW  :  Vkk  €  Pi{T),  Vr  e  Tfc,  V  €  W)  =  X  x 
=  {qe  C^(n)  :  qlK  6  Fi(K),  VK  ETh^qE  V}. 
and  the  subscript  h  for  discrete  functions  will  be  dropped  in  the  rest  of  the 
paper. 

Error  estimates  and  results  on  the  conditioning  of  the  resulting  stiffness 
matrix  can  be  found  in  [5]  (in  [15]  for  the  case  X  —  0). 

Upon  choosing  a  basis  in  and  14,  the  discrete  problem  (10)  is  turned 
into  a  linear  system  of  equations  Ax  —  b.  We  are  going  to  solve  such  system 
iteratively  using  domain  decomposition  techniques. 
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3  Domain  Decomposition  Algorithms 

We  will  introduce  and  analyze  our  domain  decomposition  algorithms  in  the 
Schwarz  framework,  which  has  been  very  successful  for  standard  Galerkin 
finite  elements,  see  [9],  [11],  [12],  [13].  We  illustrate  the  main  ideas  on  algo¬ 
rithms  that  are  representative  of  the  main  cleisses  of  domain  decomposition 
(additive,  multiplicative,  overlapping,  iterative  substructuring).  The  same 
analysis  can  be  applied  to  the  many  other  algorithms  which  have  been  pro¬ 
posed  and  analyzed  for  the  standard  scalar  case. 

We  suppose  that  the  domain  Q,  is  first  triangulated  by  a  coarse  finite 
element  triangulation  th  consisting  of  N  subdomains  0;  of  diameter  H.  The 
fine  triangulation  Th  is  a  refinement  of  th-  For  simplicity,  we  suppose  that 
each  subdomain  is  the  image  under  a  smooth  map  of  a  reference  cube. 

3.1  Overlapping  Additive  Schwarz  Methods 

Each  subdomain  Oj  is  extended  to  a  larger  subdomain  Cl'-,  consisting  of  all 
elements  of  within  a  distance  6  from  fl,  (0  <  (5  <  H). 

Each  scalar  component  of  our  finite  element  space  x  is  decomposed 
as  in  the  standard  scalar  case: 

w^  =  I:K.  =  v-‘  =  Ev:‘. 

i=l  t=l  2=1  2=1 

where 

^k,i  =  ^  ■  support(u)  C  Q'},  k  =  l,2, 3, 

Vt  =  {u^V^:  support(u)  C  fi'}. 

For  each  scalar  component,  a  global  coarse  finite  element  space  is  associated 
with  the  coarse  tri angulation  th- 

^kfi  —  ^k  =  ^  ^  trilinear  on  each  subdomain  fl,},  k  —  1,2,3, 

Vq  =  =  {p  ^  :  p  \s  trilinear  on  each  subdomain  fl,}. 

A  first  additive  method  is  defined  by  the  following  decomposition  of  the 
discrete  space,  which  maintains  the  local  and  coarse  coupling  between  the 
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different  scalar  components: 


N 

X  Wf  X 

2=0 


The  local  spaces  are 


Wj*  X  =  <  X  X 


Vy‘iXl^  i  =  l,2,--..JV 


and  the  coarse  space  is 

X  Vq  =  xV^  =  Wf  X  X  Wf  X  V^. 

We  define  the  local  projection  operators  Pi  :  x  ^  Wf  x  by 

a(Pi(u,p);  v,^)  =  a(u,p;  v,?)  V(v,q)  €  Wf  x  V/" , 
and  the  coarse  projection  operator  Po  •  x  Wq  x  Vq  by 

a(Po(u,p);v,9)  =  a(u,p;v,^)  V(v,9)  €  Wq  x  Vq. 


It  is  easy  to  see  that  the  matrix  form  of  the  local  projections  is  P  = 
RjA-^RiA,  where  the  Ri{ek)  =  |  J  otherwise'  }  '  restriction  ma¬ 

trices  selecting  only  the  unknowns  in  O'j  for  each  component  and  the  Ai  = 
RiARf  are  the  stiffness  matrices  of  local  Dirichlet  problems.  Analogously, 
Pq  =  RjjA]j  RhA,  where  Rjj  is  the  standard  piecewise  linear  interpolation 
matrix  from  the  coarse  grid  Tff  to  the  fine  grid  t/j,  for  each  component,  and 
Ah  =  RnARff  is  the  coarse  grid  discretization  of  our  problem  (9).  Let 


N 

Paddl  =  y^P- 
2=0 

The  original  discrete  problem  is  then  equivalent  to  the  preconditioned  prob¬ 
lem 

Paddlip^P^  §addl 

where  g  =  Eilo Mathew  [9].  In  matrix  form,  this 
problem  can  be  written  as  M'Mx  =  where  the  preconditioner  is 
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M~^  =  Rf  A~^Ri  +  An  optimal  convergence  bound  for  this 

algorithm  is  given  in  Theorem  2. 

A  second  additive  method  is  obtained  by  dropping  the  coupling  between 
the  different  scalar  components  of  u  and  p.  Uncoupled  local  spaces  are  now 
defined  by 

Wt-  =  X  {0}  X  {0}  X  {0}, 
w‘,.  =  {0)  X  Wl,  X  {0}  X  (0), 

Wj,.-  =  {0}  X  {0}  X  W^.  X  {0}, 

Vf  =  {0}  X  {0}  X  (0)  X  V,\ 
and  the  coarse  spaces  by 

wf  =  w‘„  =  X  {0}  X  {0}  X  {0), 

W?  =  W‘„  =  {0}  X  iy*„  X  {0}  X  {0), 

W3"  =  W3%  =  {0}  X  {0)  X  X  {0}, 

V"  =  v5  =  {0}  X  {0}  X  {0}  X 
We  then  have  the  following  decomposition 


w‘x  =  i;w‘i+f;w‘.+f;w3v+£;v'‘+wf +wf +w3«+v« 

2  =  1  2=1  2=1  2=1 


3  N 


N 


=  EEwt  +  Evf. 


As  before,  we  define  projections  Pk,i  :  x  k  =  1,2,3,  i  = 

0,  and  P4^i  :  x  U*  — >  Vf,  i  =  0,  and  the  additive 


operator 


3  AT  N 

Padd2  ~  EEa,.  +  E^4,,' 


/:=!  2=0  2=0 

We  note  that  this  algorithm  can  equivalently  be  defined  by  the  same  choice 
of  subspaces  as  for  Paddi  but  using  the  bilinear  form  b{-,-)  (introduced  in 
Theorem  1)  instead  of  a(-,  •)  in  the  definition  of  the  projections.  In  fact  this 
uncoupled  preconditioner  corresponds  to  applying  four  identical  copies  of  a 
scalar  preconditioner  to  each  scalar  component.  An  optimal  bound  also  holds 
for  this  algorithm. 


Theorem  2  There  exists  a  positive  constant  C  independent  of  h,,H  and  6 
such  that  „ 

cond(P)  <  C(1  + 

where  P  =  Paddi  or  P  =  Padd2  ■ 

Proof.  An  upper  bound  on  the  spectrum  of  P  is  standard,  since  each  point  of 
n  belongs  to  a  fixed  number  of  extended  subdomains  independent  of  N  (for 
example,  for  6  <C.  H j2  each  point  belongs  to  at  most  four  (in  2D)  or  eight  (in 
3D)  extended  subdomains).  A  lower  bound  is  obtained  by  classical  Schwarz 
analysis. 

For  P  =  Paddi,  since  we  use  exact  projections,  the  theorem  is  equivalent 
to  the  following  partition  property  (see  Dryja  and  Widlund  [13]  or  Chan  and 
Mathew  [9]): 

There  exists  a  constant  Co  such  that  V(u,p)  €  x  V'*,  there  exists  a 
decomposition  (u,p)  =  ^  ^ 

N 

^a(ui,pi;ui,pi)  <  Coo(u,p;u,p). 

izzO 

By  the  equivalence  of  Theorem  1,  this  inequality  is  equivalent  to 
N 

i=0 

which  is  a  direct  consequence  of  the  scalar  result  proven  by  Dryja  and  Wid¬ 
lund  [1-3]: 

E  wnh  <  E  ^ 

with  Cq  =  C{1  +  j). 

For  P  =  Padd2,  since  the  subspaces  are  the  same  but  we  use  inexact 
projections  defined  by  instead  of  a{-,-),  we  need  only  to  show  that 

there  exists  a  constant  u  such  that  a(u,p;u,p)  <  ujb[u,p]u,p)  V(u,p)  € 
Wf  xVf-,  i  =  0, 1,---,A  (see  Dryja  and  Widlund  [12]).  This  follows 
immediately  from  the  equivalence  of  a  and  b. 
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3.2  Overlapping  Multiplicative  Schwarz  Methods 

By  using  the  same  coupled  local  and  coarse  spaces  as  in  the  additive  algorithm 
Padd\  1  we  can  define  a  multiplicative  operator: 

PmuH  =  I  -  {I  -  -  Px){I  -  Po). 

The  multiplicative  algorithm  consists  in  solving  the  nonsymmetric  system 


Pmult{}^iP)  —  Sm,ult 

by  an  iterative  method  such  as  GMRES. 

We  can  also  define  a  symmetrized  multiplicative  operator 


PmuHs  =  /-(/-  Po)  •••(/  -  Pn-i){I  -  Pn){I  -  P;V-i)  •  •  •  (/  -  Po) 


and  a  symmetrized  algorithm,  consisting  in  solving  the  symmetric  system 


P multsiy^i  p}  —  ^mults 
by  an  iterative  method  like  CG. 

Theorem  3  There  exists  a  positive  constant  C  independent  of  h,  H  and  8 
such  that  ^ 

COnd{Pmults)  <  ^^(1  +  y). 

The  proof  is  again  based  on  the  extension  of  the  scalar  result  (see,  for  exam¬ 
ple,  Chan  and  Mathew  [9])  by  using  the  equivalence  of  Theorem  1.  Analo¬ 
gously,  multiplicative  versions  of  Padd2  could  be  built  using  uncoupled  local 
and  coarse  spaces. 


3.3  An  Iterative  Substructuring  Method 

For  a  complete  and  detailed  analysis  of  this  class  of  methods,  we  refer  to 
Dryja,  Smith  and  Widlund  [11].  Here  we  consider  only  a  simple  represen¬ 
tative  of  this  class,  namely  the  analog  of  Algorithm  6.2  in  [11],  which  is 
vertex-based  and  has  a  standard  coarse  space.  For  simplicity,  we  only  con¬ 
sider  the  uncoupled  additive  version. 
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The  standard  first  step  of  nonoverlapping  methods  is  the  elimination  of 
the  variables  interior  to  each  subdomain  (at  least  implicitly).  We  then  work 
with  the  Schur  complement  S  =  Kbb  —  Kib  of  the  stiffness  matrix 


K  = 


Kn  Kib  \ 
KJb  Kbb  j  ' 


The  reduced  linear  system  with  S  involves  only  variables  on  the  interface 
r  =  \  Ti).  When  solving  with  a  preconditioned  iterative  method,  we 

need  only  the  action  of  5"  on  a  given  vector  and  there  is  no  need  to  assemble 
S  explicitly. 

In  the  Schwarz  framework,  working  with  S  corresponds  to  working  with 
the  discrete  harmonic  subspace  x  of  the  original  space  x  V^. 
Local  spaces  are  associated  with  the  geometric  objects  (faces  edges  Ei 
and  vertices  Vi)  forming  the  interface  F.  Each  scalar  space  is  decomposed  as 

Wt  =  E  W'.V,  +  E  wt  E,  +  E  H.. .  *  =  1.2.3. 

Fi  Ei  Vi 


and 

v^=i:ni+Eyk+En- 

Fi  Ei  vi 

Here,  for  example,  W^p.  =  {u  G  :  u  =  0  on  Th  —  Fi^h}-,  where  F/i  and  Fi^h 
are  the  set  of  nodes  on  F  and  Fi  respectively.  The  other  spaces  are  defined 
analogously.  As  for  the  overlapping  case,  we  then  embed  these  scalar  spaces 
in  our  product  space  x  for  example,  =  W^^p.  x{0}x{0}x{0}. 
As  a  coarse  space,  we  consider  the  discrete  harmonic  subspace  of  the  same 
coarse  space  used  for  Padd2,  i-e.,  +  W|^  +  V^.  We  obtain  the 

following  decomposition 

W*  X  i/‘  =  E(Ew{.f,  +  E  +  Ewt,  +  wf) 

k=l  Fi  Ei  Vi 


+ev^.,+ev6+ev;,+v". 

Fi  Ei  Vi 

By  defining  as  before  projection  operators  into  the  subspaces,  we  form  the 
additive  operator 

Pis  =  +  X]  -^^•3'.  +  ^k,o), 

k=l  Fi  Ei  Vi 
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where  again  for  /;  =  4  the  projections  are  into  the  Vf  spaces. 

Theorem  4  There  exists  a  positive  constant  C  independent  of  h  and  H  such 
that 

cond{Pis)  <  C{1  +  \og{H/h)Y. 

As  before,  the  proof  is  based  on  the  extension  of  the  scalar  result  (see  Dryja, 
Smith  and  Widlund  [11],  Theorem  6.2)  by  using  the  equivalence  of  Theorem 
1. 

4  Numerical  Results 

In  this  section,  we  report  the  results  of  numerical  experiments  which  confirm 
the  optimal  convergence  bounds  obtained  in  the  previous  sections.  All  the 
results  have  been  obtained  with  Matlab  4.2  running  on  Sun  Sparcstations. 
The  model  problem  considered  is  the  standard  Poisson  equation  {A  —  I,  X  = 
0)  on  the  unit  square,  with  p  =  0  on  Pd  =  and  77.U  =  0  on  (i.e.  ui  =  0 
on  {y  =  0}  and  [y  =  1);  U2  =  0  on  {a;  =  0}  and  [x  =  1}).  The  right-hand 
side  /  is  chosen  such  that  we  have  p(x,  y)  =  sin{'Kx)sin{'Ky)  as  exact  solution. 
0  is  decomposed  into  a  regular  grid  of  N  square  subdomains,  with  N  varying 
from  2  X  2  to  8  X  8.  The  fine  grid  mesh  size  h  varies  from  1/.32  to  1/128. 

The  Krylov  method  used  for  all  the  symmetric  problems  is  PCG,  while 
we  use  GMRES  for  the  nonsymmetric  problem  with  Pmuit  ■  The  initial  guess 
is  always  zero  and  the  stopping  criterion  is  ||rfc||2/||ro||2  <  10“®,  where  Vk  is 
the  residual  at  step  h. 

The  local  and  coarse  problems  involved  in  the  application  of  the  precon¬ 
ditioners  are  always  solved  directly.  For  each  method,  we  report  the  number 
of  iterations  and  Lanczos-based  estimates  of  the  condition  number  and  the 
extreme  eigenvalues  (except  for  the  multiplicative  algorithm,  where  we  report 
the  average  convergence  factor  instead). 

Overlapping  additive  methods.  We  have  first  studied  the  coupled  method 
Paddi  with  fixed  minimal  overlap  size  8  =  h.  The  mesh  size  h  is  decreased 
while  the  number  of  subdomains  N  is  increased  proportionally,  so  that  the 
subdomain  size  H/h  =  is  kept  constant  {H  =  l/\/]V).  The  results  are  re¬ 
ported  in  Table  1  and  clearly  show  a  constant  condition  number  cond{Paddi )  = 
Xmaxl^mim  for  problem  sizes  from  3007(iV  =  4)  to  48895(iV  =  64). 
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Table  1:  Paddi'-  Overlapping  Additive  Schwarz  with  fixed  overlap  size  6  =  h. 


N 

h-^ 

iter. 

cond{Padd\) 

^max 

^min 

4 

32 

16 

11.2172 

4.0048 

0.3570 

9 

48 

19 

12.1787 

4.0068 

0.3290 

16 

64 

20 

11.9775 

4.0050 

0.3343 

25 

80 

20 

11.1689 

4.0052 

0.3586 

36 

96 

21 

12.5450 

4.0044 

0.3192 . 

49 

112 

20 

11.9944 

4.0050 

0.3339 

64 

128 

21 

12.5500 

4.0047 

0.3191 

Table  2:  Paddi-  Overlapping  Additive  Schwarz  with  fixed  number  of  subdo¬ 
mains  N  =  6T _ _ _ _ 


6 

h-^ 

iter. 

cond{Paddi) 

^max 

^min 

h 

128 

21 

12.5500 

4.0047 

0.3191 

2h 

128 

17 

7.1316 

4.0307 

0.5651 

3h 

128 

16 

5.5769 

4.0765 

0.7309 

4h 

128 

15 

4.9540 

4.1396 

0.8356 

5h 

128 

15 

4.6460 

4.2170 

0.9076 

6h 

128 

15 

4.5125 

4.3054 

0.9541 

7h 

128 

16 

4.5859 

4.4018 

0.9598 
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Table  3:  Padd2  ’■  Overlapping  Additive  Schwarz  with  fixed  overlap  size  6  =  h. 


N 

iter. 

cond{Padd2) 

^max 

^min 

4 

32 

17 

10.3521 

4.0050 

0.3868 

9 

48 

20 

12.6290 

4.0051 

0.3171 

16 

64 

20 

11.9811 

4.0051 

0.3342 

25 

80 

21 

11.3821 

4.0052 

0.3518 

36 

96 

21 

12.5458 

4.0043 

0..3191 

49 

112 

20 

11.9997 

4.0052 

0.3337 

64 

128 

21 

12.5261 

4.0047 

0.3197 

Table  4:  Padd2  '■  Overlapping  Additive  Schwarz  with  fixed  number  of  subdo¬ 
mains  N  =  64.  _ _ 


6 

h-^ 

iter. 

cond{Padd2) 

^max 

^min 

h 

128 

21 

12.5261 

4.0047 

0.3197 

2h 

128 

17 

7.1206 

4.0315 

0.5661 

3h 

128 

16 

5.5513 

4.0777 

0.7345 

4h 

128 

16 

5.3850 

4.1442 

0.7695 

5h 

128 

16 

5.4545 

4.2233 

0.7742 

6h 

128 

16 

5.5306 

4.3158 

0.7803 

7h 

128 

16 

5.6176 

4.4297 

0.7885 
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Table  5:  Pmuit  and  Pmults'-  Overlapping  Multiplicative  Schwarz  with  fixed 
overlap  size  8  =  h. 


N 

h-^ 

multiplicative  (GMRES) 
iter.  p={rilroY^^ 

symmetrized  multiplicative  (CG) 
iter.  COnd(^Pmults^  ^max  ^min 

4 

32 

8 

0.1847 

7 

1.8576 

0.9994 

0.5379 

9 

48 

7 

0.1433 

6 

1.7398 

0.9999 

0.5749 

16 

64 

6 

0.1233 

6 

1.7600 

0.9999 

0.5681 

25 

80 

6 

0.1102 

6 

1.6810 

0.9999 

0.5948 

36 

96 

6 

0.1021 

6 

1.6940 

0.9999 

0.5902 

49 

112 

6 

0.0952 

6 

1.6661 

0.9999 

0.6001 

64 

128 

5 

0.0849 

6 

1.7308 

0.9999 

0.6079 

In  Table  2,  we  fix  the  mesh  size  [h  =  1/128)  and  the  decomposition 
(iV  =  64)  and  we  vary  the  overlap  size  8  from  h  to  7h.  As  in  the  scalar  case, 
the  condition  number  cond{Paddi)  improves  as  8  increases,  because  of  X^in 
being  closer  to  unity.  For  large  overlap,  the  improvement  becomes  negligible 
or  negative,  because  of  the  growth  of  Xmax- 

The  same  sets  of  results  for  the  uncoupled  method  Padd2  are  reported 
in  Table  3  and  Table  4,  respectively.  For  this  simple  model  problem,  the 
uncoupled  method  is  only  slightly  worse  than  the  coupled  one,  in  terms  of 
iteration  count  (some  condition  numbers  are  almost  the  same  or  even  better 
for  Padd2)-  We  point  out  that  although  A  =  /,  eliminating  diffusive  coupling 
between  the  flux  components,  there  is  still  coupling  between  the  flux  variables 
and  p,  so  the  strong  performance  of  Padd2  is  encouraging. 

Overlapping  multiplicative  methods.  In  Table  5,  we  compare  the  multi¬ 
plicative  method  Pmuit  accelerated  with  GMRES  and  the  symmetrized  mul¬ 
tiplicative  method  Pmults  accelerated  with  CG.  We  consider  the  two  methods 
with  minimal  overlap  and  constant  subdomain  size.  Since  Pmuit  is  nonsym- 
metric,  we  report  the  average  convergence  factor  p  =  (rj/ro)^^*  instead  of  the 
condition  number.  Even  if  the  symmetrized  version  is  approximately  twice 
as  expensive  as  the  standard  one,  the  number  of  iterations  is  almost  the  same 
for  the  two  methods.  Therefore,  the  symmetrized  version  is  less  efficient  on 
this  simple  problem. 
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Table  6:  Pis'.  Iterative  Substructuring  . 


N 

h-^ 

iter. 

cond{Pis) 

^max 

^min 

4 

32 

9 

3.40.35 

1.5691 

0.4610 

9 

48 

17 

7.8812 

1.8497 

0.2347 

16 

64 

18 

7.8543 

1.7962 

0.2287 

25 

80 

18 

8.5822 

1.8864 

0.2198 

36 

96 

19 

9.4115 

1.8511 

0.1966 

49 

112 

18 

8.6646 

1.8939 

0.2185 

64 

128 

19 

9.6532 

1.8617 

0.1928 

Iterative  substructuring.  Table  6  shows  the  results  for  the  iterative  sub- 
structuring  methods  Pis  with  fixed  subdomain  size.  They  clearly  show  a 
constant  bound  for  the  condition  number  and  the  number  of  iterations. 

5  Conclusions 

In  this  paper,  some  domain  decomposition  algorithms  have  been  introduced 
for  the  discrete  systems  arising  from  first-order  system  least  squares  methods 
applied  to  second-order  elliptic  problems.  These  recently  proposed  methods 
allow  the  use  of  standard  finite  element  spaces,  which  are  not  required  to 
satisfy  the  inf-sup  condition. 

The  analysis  of  the  domain  decomposition  algorithms  follows  from  analo¬ 
gous  results  for  the  standard  Galerkin  case  and  the  equivalence  between  the 
bilinear  form  associated  with  the  least  squares  functional  and  the 
norm. 

Optimal  convergence  bounds  have  been  proven  for  overlapping  algorithms 
(additive,  multiplicative,  coupled,  uncoupled  versions),  while  quasi-optimal 
bounds  have  been  proven  for  iterative  substructuring  algorithms.  Numerical 
experiments  on  a  simple  model  problem  confirm  these  bounds. 

Future  work  will  investigate  the  performance  of  these  algorithms  for  prob¬ 
lems  with  convection  and  for  elliptic  systems. 
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